This document present some examples of vector field processing using pygsf methods.
In order to plot fields, we run the following commands:
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
The plot codes are modified from [1] as answered by nicoguaro.
The modules to import for dealing with grids are:
In [2]:
from pygsf.mathematics.arrays import *
from pygsf.spatial.rasters.geotransform import *
from pygsf.spatial.rasters.fields import *
The definition of divergence for our 2D case is:
The definition of curl module in our 2D case is:
so that the module of the curl is:
The implementation of the curl module calculation has been debugged using the code at [2] by Johnny Lin. Deviations from the expected theoretical values are the same for both implementations.
We calculate a theoretical, 2D vector field and check that the parameters calculated by pygsf is equal to the expected one.
We use a modified example from p. 67 in [3].
In order to create the two grids that represent the x- and the y-components, we therefore define the following two "transfer" functions from coordinates to z values:
In [3]:
def z_func_fx(x, y):
return 0.0001 * x * y**3
def z_func_fy(x, y):
return - 0.0002 * x**2 * y
The above functions define the value of the cells, using the given x and y geographic coordinates.
Gridded field values are calculated for the theoretical source vector field x- and y- components using the provided number of rows and columns for the grid:
In [4]:
rows=100; cols=200
In [5]:
size_x = 5; size_y = 5
In [6]:
tlx = 500.0; tly = 250.0
Arrays components are defined in terms of indices i and j, so to transform array indices to geographical coordinates we use a geotransform. The one chosen is:
In [7]:
gt1 = GeoTransform(
inTopLeftX=tlx,
inTopLeftY=tly,
inPixWidth=size_x,
inPixHeight=size_y)
Note that the chosen geotransform has no axis rotation, as is in the most part of cases with geographic grids.
In [8]:
fx1 = array_from_function(
row_num=rows,
col_num=cols,
geotransform=gt1,
z_transfer_func=z_func_fx)
In [9]:
print(fx1)
In [10]:
fy1 = array_from_function(
row_num=rows,
col_num=cols,
geotransform=gt1,
z_transfer_func=z_func_fy)
In [11]:
print(fy1)
To visualize the parameters of the flow, we calculate the geographic coordinates:
In [12]:
X, Y = gtToxyCellCenters(
gt=gt1,
num_rows=rows,
num_cols=cols)
and the vector field magnitude:
In [13]:
magn = magnitude(
fld_x=fx1,
fld_y=fy1)
We can then visualize it:
In [14]:
plt.figure(figsize=(14, 6))
plt.contourf(X, Y, np.log10(magn), cmap="bwr")
cbar = plt.colorbar()
cbar.ax.set_ylabel('Magnitude (log10)')
plt.streamplot(X, Y, fx1/magn, fy1/magn, color="black")
plt.axis("image")
plt.title('Vector field magnitude and streamlines')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
Since the vector field formula is:
the theoretical divergence transfer function is:
In [15]:
def z_func_div(x, y):
return 0.0001 * y**3 - 0.0002 * x**2
The theoretical divergence field can be created using the function expressing the analytical derivatives z_func_div:
In [16]:
theor_div = array_from_function(
row_num=rows,
col_num=cols,
geotransform=gt1,
z_transfer_func=z_func_div)
In [17]:
print(theor_div)
In [18]:
plt.figure(figsize=(14, 6))
plt.contourf(X, Y, theor_div, cmap="RdYlGn_r")
cbar = plt.colorbar()
cbar.ax.set_ylabel('Theoretical divergence')
plt.streamplot(X, Y, fx1/magn, fy1/magn, color="black")
plt.axis("image")
plt.title('Vector field theoretical divergence and streamlines')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
Divergence as resulting from pygsf calculation:
In [19]:
div = divergence(
fld_x=fx1,
fld_y=fy1,
cell_size_x=size_x,
cell_size_y=size_y)
In [20]:
print(div)
In [21]:
plt.figure(figsize=(14, 6))
plt.contourf(X, Y, div, cmap="RdYlGn_r")
cbar = plt.colorbar()
cbar.ax.set_ylabel('Empirical divergence')
plt.streamplot(X, Y, fx1/magn, fy1/magn, color="black")
plt.axis("image")
plt.title('Vector field empirical divergence and streamlines')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
We check whether the theoretical and the estimated divergence fields are close:
In [22]:
np.allclose(theor_div, div)
Out[22]:
The vector function is:
therefore the theoretical curl module is:
so that the theoretical transfer function is:
In [23]:
def z_func_curl_mod(x, y):
return - 0.0004 * x * y - 0.0003 * x * y**2
The theoretical divergence field can be created using the function expressing the analytical derivatives z_func_div:
In [24]:
theor_curl_mod = array_from_function(
row_num=rows,
col_num=cols,
geotransform=gt1,
z_transfer_func=z_func_curl_mod)
In [25]:
print(theor_curl_mod)
In [26]:
plt.figure(figsize=(14, 6))
plt.contourf(X, Y, theor_curl_mod, cmap="spring")
cbar = plt.colorbar()
cbar.ax.set_ylabel('Theoretical curl module')
plt.streamplot(X, Y, fx1/magn, fy1/magn, color="black")
plt.axis("image")
plt.title('Vector field theoretical curl module and streamlines')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
The module of curl as resulting from pygsf calculation is:
In [27]:
curl_mod = curl_module(
fld_x=fx1,
fld_y=fy1,
cell_size_x=size_x,
cell_size_y=size_y)
In [28]:
print(curl_mod)
In [29]:
plt.figure(figsize=(14, 6))
plt.contourf(X, Y, theor_curl_mod, cmap="spring")
cbar = plt.colorbar()
cbar.ax.set_ylabel('Empirical curl module')
plt.streamplot(X, Y, fx1/magn, fy1/magn, color="black")
plt.axis("image")
plt.title('Vector field empirical curl module and streamlines')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
We check whether the theoretical and the estimated curl module fields are close:
In [30]:
np.allclose(theor_curl_mod, curl_mod)
Out[30]:
We look at where there are significant differences between the theoretical and the empiric curl module fields, by calculating the percent difference between these two fields:
In [31]:
percent_diffs = 100.0*(curl_mod - theor_curl_mod)/theor_curl_mod
plt.figure(figsize=(14, 6))
plt.contourf(X, Y, percent_diffs, cmap="gist_stern_r")
cbar = plt.colorbar()
cbar.ax.set_ylabel('Difference empirical-theoretical curl module')
plt.streamplot(X, Y, fx1/magn, fy1/magn, color="black")
plt.axis("image")
plt.title('Differences between empirical and theoretical curl module')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
While the majority of the empiric values are near the theoretical ones, there is an almost horizontal strip at around y = 0 where flows have singularities, with very high deviances (both negative and positive) that are related to the observed "singularity" of the v_x field gradient along the y axis.
In [32]:
plt.hist(percent_diffs) # arguments are passed to np.histogram
plt.show()
We test another theoretical, 2D vector field, maintaining the same geotransform and other grid parameters as in the previous example. We use the field described in example 1 in [4]:
The "transfer" functions from coordinates to z values are:
In [33]:
def z_func_fx(x, y):
return y
def z_func_fy(x, y):
return - x
Gridded field values are calculated for the theoretical source vector field x- and y- components using the provided number of rows and columns for the grid:
In [34]:
rows=200; cols=200
In [35]:
size_x = 10; size_y = 10
In [36]:
tlx = -1000.0; tly = 1000.0
Arrays components are defined in terms of indices i and j, so to transform array indices to geographical coordinates we use a geotransform. The one chosen is:
In [37]:
gt1 = GeoTransform(
inTopLeftX=tlx,
inTopLeftY=tly,
inPixWidth=size_x,
inPixHeight=size_y)
Note that the chosen geotransform has no axis rotation, as is in the most part of cases with geographic grids.
In [38]:
fx2 = array_from_function(
row_num=rows,
col_num=cols,
geotransform=gt1,
z_transfer_func=z_func_fx)
In [39]:
print(fx2)
In [40]:
fy2 = array_from_function(
row_num=rows,
col_num=cols,
geotransform=gt1,
z_transfer_func=z_func_fy)
In [41]:
print(fy2)
We visualize the parameters of the flow.
The geographic coordinates are:
In [42]:
X, Y = gtToxyCellCenters(
gt=gt1,
num_rows=rows,
num_cols=cols)
The vector field magnitude is:
In [43]:
magn = magnitude(
fld_x=fx2,
fld_y=fy2)
In [44]:
plt.figure(figsize=(14, 14))
plt.contourf(X, Y, magn, cmap="cool")
plt.colorbar()
plt.streamplot(X, Y, fx2/magn, fy2/magn, color="black")
plt.axis("image")
plt.show()
The theoretical curl module is a constant value:
The module of curl as resulting from pygsf calculation is:
In [45]:
curl_mod = curl_module(
fld_x=fx2,
fld_y=fy2,
cell_size_x=size_x,
cell_size_y=size_y)
In [46]:
print(curl_mod)
In [47]:
plt.figure(figsize=(14, 10))
plt.contourf(X, Y, curl_mod, cmap="bwr")
plt.colorbar()
plt.streamplot(X, Y, fx2/magn, fy2/magn, color="black")
plt.axis("image")
plt.show()
We check whether the theoretical and the estimated curl module fields are close:
In [48]:
np.allclose(-2.0, curl_mod)
Out[48]:
Deviances from the expected value are reduced, as evidenced by the histogram:
In [49]:
plt.hist(curl_mod + 2.0) # arguments are passed to np.histogram
plt.show()
[1] Visually appealing ways to plot singular vector fields with matplotlib or other foss tools. https://scicomp.stackexchange.com/questions/18760/visually-appealing-ways-to-plot-singular-vector-fields-with-matplotlib-or-other
[2] Curl Function (Solution). http://www.johnny-lin.com/ams2011/sc/arrays_io/as/hide/py-curl-soln.shtml. Consulted on June 4, 2018.
[3] M. R. Spiegel, 1975. Analisi Vettoriale. Etas Libri, pp. 224.
[4] Curl (mathematics). https://en.wikipedia.org/wiki/Curl_(mathematics). Consulted on June 6, 2018.